setwd("")
library(foreign)
library(dplyr)
library(mgcv)
source("functions.R")

################################## ################################## 
################################## ################################## 
############ Appendix Figures
################################## ################################## 
################################## ################################## 


################################## ################################## 
################# Figure B1 - GAM

pdf("Figures/FigureB1.pdf",width=8, height=4)
par(mfrow=c(1,2))

swissdata <- read.dta("./Datasets/Other/Swiss_master.dta")
swissdata <- subset(swissdata,(reformyear<=0))
swissdata <- subset(swissdata, pop > quantile(swissdata$pop,.05,na.rm=T) & pop < quantile(swissdata$pop,.95,na.rm=T)) # Population trim
swissdata$interaction_constant <- swissdata$interaction_constant*100
swissdata$log_net_welfare_head <- swissdata$log_net_welfare_head*100

mod.gam <- gam(log_net_welfare_head ~  s(interaction_constant,sp=5) +year + regunemployed_pop + lowbracket_rate + sp + log_rev + logpop + pctforeign + factor(bfsnr), data=swissdata)
plot(mod.gam,se=1.64, main="Switzerland", xlab="% Foreign Residents",shade=T,cex.lab=.75,cex.axis=.7,shade.col="gray90",ylab="% Change in Social Spending",
     rug=T,xlim=c(0,quantile(swissdata$interaction_constant[swissdata$reformyear==0],c(.975),na.rm=T)))
grid()

bdata<- read.dta("./Datasets/Other/Belgium_master.dta")
bdata <- subset(bdata,bdata$year <= 2005)
bdata$pctnoneunoncit <- bdata$pctnoneunoncit*100
bdata$main_interaction <- bdata$main_interaction*100
bdata$log_cpas_head <- bdata$log_cpas_head*100
mod.gam <- gam(log_cpas_head~  s(main_interaction,sp=2) + factor(year) + pctnoneunoncit + logpop + filing_rate + unemp_rate + socshare + log_rev + factor(id), data=bdata)
plot(mod.gam,se=1.64, main="Belgium", xlab="% Non-EU Foreign Residents",cex.lab=.75,cex.axis=.7,shade=T,shade.col="gray90",ylab="% Change in Social Spending",
     xlim=c(0,quantile(bdata$pctnoneunoncit,c(.975),na.rm=T)))
grid()
dev.off()

################################## ################################## 
################# Figure B2 - LOESS

pdf("Figures/FigureB2.pdf",width=8, height=4)
par(mfrow=c(1,2))

data<- read.dta("./Datasets/Other/Swiss_master.dta")
data <- subset(data, data$reformyear==2) 
data <- subset(data, data$pop > quantile(data$pop,.05,na.rm=T) & pop < quantile(data$pop,.95,na.rm=T)) # Remove population outiers
data <- subset(data, I4_net_welfare_head < quantile(data$I4_net_welfare_head ,c(.975),na.rm=T) & I4_net_welfare_head > quantile(data$I4_net_welfare_head ,c(.025),na.rm=T)) 
x <- data$pctforeign*100
y <- data$I4_net_welfare_head*100
general_loess(x,y,1,"% Foreign Residents","Change in Welfare Spending (Indexed to 100)",c(120,180),"Switzerland",0,.975)

belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- subset(belgiumdata,year==2007)
belgiumdata <- subset(belgiumdata, I4_cpas_head < quantile(belgiumdata$I4_cpas_head ,c(.975),na.rm=T) & I4_cpas_head > quantile(belgiumdata$I4_cpas_head ,c(.025),na.rm=T)) 
x <- belgiumdata$pctnoneunoncit*100
y <- belgiumdata$I4_cpas_head*100
general_loess(x,y,1,"% Non-EU Foreign Residents","Change in Welfare Spending (Indexed to 100)",c(100,130),"Belgium",0,.99)
dev.off()

################################## ################################## 
################# Figure B3 - Unemployment

pdf("Figures/FigureB3.pdf",width=8, height=4)
par(mfrow=c(1,2))

############ Switzerland
swissdata <- read.dta("./Datasets/Other/Swiss_master.dta")
data <- swissdata %>% arrange(bfsnr,year)  %>% group_by(bfsnr) %>% mutate(d1_unemp = (regunemployed_pop - lag(regunemployed_pop,1)) ) 
data <- subset(data,data$diff_voteright_following == 1)
x <- (data$pctforeign)*100
y <- data$d1_unemp*100
general_loess(x,y,1,"% Foreign Residents","Change in Unemployment Rate (Percentage Points)",c(-.5,1),"Switzerland",0,.975)
abline(h=0,col="gray")


############ Belgium
belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- belgiumdata %>% arrange(id,year)  %>% group_by(id) %>% mutate(d1_unemp = (unemp_rate - lag(unemp_rate,1)) ) 
belgiumdata <- subset(belgiumdata,belgiumdata$year==2005)
belgiumdata <- subset(belgiumdata, pctnoneunoncit < quantile(belgiumdata$pctnoneunoncit, c(.975),na.rm=T))
data <- belgiumdata
x <- (data$pctnoneunoncit)*100
y <- data$d1_unemp
general_loess(x,y,.9,"% Foreign Residents","Change in Unemployment Rate (Percentage Points)",c(-.5,1),"Belgium",0,.975)
abline(h=0,col="gray")
dev.off()

################################## ################################## 
################# Figure B4 - RIS Caseload

pdf("Figures/FigureB4.pdf",width=6, height=4)
belgiumdata<- read.dta("./Datasets/Other/Belgium_master.dta")
belgiumdata <- belgiumdata %>% arrange(id,year)  %>% group_by(id) %>% mutate(l1_ris_enrollment = lag(ris_enrollment,1)) 
belgiumdata <- subset(belgiumdata,belgiumdata$year==2005)
belgiumdata <- subset(belgiumdata, pctnoneunoncit < quantile(belgiumdata$pctnoneunoncit, c(.99),na.rm=T))
belgiumdata$ris_change <- (belgiumdata$ris_enrollment - belgiumdata$l1_ris_enrollment)/belgiumdata$l1_ris_enrollment
belgiumdata$ris_change[!is.finite(belgiumdata$ris_change)] <- NA
mdata <- na.omit(data.frame(belgiumdata$pctnoneunoncit,belgiumdata$ris_change))
x <- (mdata[,1])*100
y <- (mdata[,2])*100
general_loess(x,y,.9,"% Non-EU","% Change in RIS Caseload",c(-25,30),"")
abline(h=0,col="gray")
dev.off()

################################## ################################## 
################# Figure B5 - GAM Individual Cantons
pdf("Figures/FigureB5.pdf",width=8, height=4)
par(mfrow=c(1,3))

swissdata <- read.dta("./Datasets/Other/Swiss_master.dta")
swissdata <- subset(swissdata,(reformyear<=0))
swissdata <- subset(swissdata, pop > quantile(swissdata$pop,.05,na.rm=T) & pop < quantile(swissdata$pop,.95,na.rm=T)) # Population trim
swissdata$interaction_constant <- swissdata$interaction_constant*100
swissdata$log_net_welfare_head <- swissdata$log_net_welfare_head*100

swissdata2 <- subset(swissdata, cantonid==10) 
mod.gam <- gam(log_net_welfare_head ~  s(interaction_constant,sp=1) + year + regunemployed_pop + lowbracket_rate + sp + log_rev + logpop + pctforeign + factor(bfsnr), data=swissdata2)
plot(mod.gam,se=1.64,main="Fribourg", xlab="% Foreign Residents",shade=T,cex.lab=.75,cex.axis=.7,shade.col="gray95",ylab="% Change in Social Spending",
     rug=T,xlim=c(100*min(swissdata2$pctforeign),quantile(swissdata2$interaction_constant[swissdata2$reformyear==0],c(.95),na.rm=T)))
grid()

swissdata2 <- subset(swissdata, cantonid==25) 
mod.gam <- gam(log_net_welfare_head ~  s(interaction_constant,sp=.15) + year + regunemployed_pop + lowbracket_rate + sp + log_rev + logpop + pctforeign + factor(bfsnr), data=swissdata2)
plot(mod.gam,se=1.64,main="Geneva", xlab="% Foreign Residents",shade=T,cex.lab=.75,cex.axis=.7,shade.col="gray95",ylab="% Change in Social Spending",
     rug=T,xlim=c(100*min(swissdata2$pctforeign),quantile(swissdata2$interaction_constant[swissdata2$reformyear==0],c(.95),na.rm=T)))
grid()

swissdata2 <- subset(swissdata, cantonid==22) 
mod.gam <- gam(log_net_welfare_head ~  s(interaction_constant,sp=3) + year + regunemployed_pop + lowbracket_rate + sp + log_rev + logpop + pctforeign + factor(bfsnr), data=swissdata2)
plot(mod.gam,se=1.64, main="Vaud", xlab="% Foreign Residents",shade=T,cex.lab=.75,cex.axis=.7,shade.col="gray95",ylab="% Change in Social Spending",
     rug=T,xlim=c(100*min(swissdata2$pctforeign),quantile(swissdata2$interaction_constant[swissdata2$reformyear==0],c(.95),na.rm=T)))
grid()
dev.off()


################################## ################################## 
################# Figure B6 - LOESS Individual Cantons

pdf("Figures/FigureB6.pdf",width=8, height=4)
par(mfrow=c(1,3))
data<- read.dta("./Datasets/Other/Swiss_master.dta")
data <- subset(data, data$reformyear==2) 
data <- subset(data, data$pop > quantile(data$pop,.05,na.rm=T) & pop < quantile(data$pop,.95,na.rm=T)) # Population outliers

dataf <- subset(data, cantonid==10) 
x <- dataf$pctforeign*100
y <- (dataf$s4_net_welfare_head )
general_loess(x,y,1,"% Foreign Residents","Change in Welfare Spending (Indexed to 100)",c(-10,200),"Fribourg",.025,.975)

datag <- subset(data, cantonid==25) 
datag <- subset(datag, datag$pctforeign > quantile(datag$pctforeign,c(.025),na.rm=T) & datag$pctforeign < quantile(datag$pctforeign,c(.975),na.rm=T)) 
x <- datag$pctforeign*100
y <- (datag$s4_net_welfare_head)
general_loess(x,y,1,"% Foreign Residents","Change in Welfare Spending (Indexed to 100)",c(-10,175),"Geneva",.025,.975)

datav <- subset(data, cantonid==22) 
datav <- subset(datav, datav$pctforeign > quantile(datav$pctforeign,c(.025),na.rm=T) & datav$pctforeign < quantile(datav$pctforeign,c(.975),na.rm=T)) 
x <- datav$pctforeign*100
y <- (datav$s4_net_welfare_head )
general_loess(x,y,1,"% Foreign Residents","Change in Welfare Spending (Indexed to 100)",c(0,300),"Vaud",.025,.975)

dev.off()


################################## ################################## 
################# Figure B7 - Revenue

pdf("Figures/FigureB7.pdf",width=8, height=4)
par(mfrow=c(1,2))

########### SWITZERLAND

sdata <- read.dta("./Datasets/Other/Swiss_master.dta")
sdata <- subset(sdata,reformyear==2)  
sdata <- subset(sdata,rev_3year > quantile(sdata$rev_3year,.025,na.rm=T) & rev_3year < quantile(sdata$rev_3year,.975,na.rm=T))
mdata <- na.omit(data.frame(sdata$rev_3year,sdata$interaction_constant))
x <- mdata$sdata.interaction_constant*100
y <- mdata$sdata.rev_3year - mean(mdata$sdata.rev_3year,na.rm=T)
y <- y*100
general_loess(x,y,1,"% Foreign Residents in Year of Reform","% Change in Local Revenue",c(-15,25),"Switzerland",0,.99)

# Placebo
sdata <- read.dta("./Datasets/Other/Swiss_master.dta")
sdata <- subset(sdata,reformyear==-1)  
sdata <- subset(sdata,rev_3year > quantile(sdata$rev_3year,.025,na.rm=T) & rev_3year < quantile(sdata$rev_3year,.975,na.rm=T))
sdata <- subset(sdata, pctforeign < quantile(sdata$pctforeign,.99,na.rm=T)) 
mdata <- na.omit(data.frame(sdata$rev_3year,sdata$pctforeign_fixed))
x <- mdata$sdata.pctforeign_fixed*100
y <- mdata$sdata.rev_3year - mean(mdata$sdata.rev_3year,na.rm=T)
y <- y*100
lo <- loess(y~x,span=1)
xl <- seq(min(x),max(x), (max(x)- min(x))/1000)
pred.c <- predict(lo,xl,se=T)
lines(xl, pred.c$fit, col='black', lty=2)
legend("topleft", c("2 Years Prior to Reform", "2 Years Following Reform"),lty=c(2,1),lwd=c(1,1),cex=.6,bg="white")


########### BELGIUM

bdata<- read.dta("./Datasets/Other/Belgium_master.dta")
bdata <- subset(bdata,year==2007)  
bdata <- subset(bdata,rev_3year > quantile(bdata$rev_3year,.05,na.rm=T) & rev_3year < quantile(bdata$rev_3year,.95,na.rm=T))
mdata <- na.omit(data.frame(bdata$rev_3year,bdata$main_interaction))
x <- mdata$bdata.main_interaction*100
y <- mdata$bdata.rev_3year - mean(mdata$bdata.rev_3year,na.rm=T)
y <- y*100
general_loess(x,y,1,"% Non-EU Residents in Year of Reform","% Change in Local Revenue",c(-7.5,10),"Belgium",0,.99)

# Placebo 
bdata<- read.dta("./Datasets/Other/Belgium_master.dta")
bdata <- subset(bdata,year==2004)  
bdata <- subset(bdata,pctnoneunoncit < quantile(bdata$pctnoneunoncit,.99,na.rm=T))
bdata <- subset(bdata,rev_3year > quantile(bdata$rev_3year,.05,na.rm=T) & rev_3year < quantile(bdata$rev_3year,.95,na.rm=T))
mdata <- na.omit(data.frame(bdata$rev_3year,bdata$interaction_fixed))
x <- mdata$bdata.interaction_fixed*100
y <- mdata$bdata.rev_3year - mean(mdata$bdata.rev_3year,na.rm=T)
y <- y*100
lo <- loess(y~x,span=1)
xl <- seq(min(x),max(x), (max(x) - min(x))/1000)
pred.c <- predict(lo,xl,se=T)
lines(xl, pred.c$fit, col='black', lty=2)
legend("topleft", c("2 Years Prior to Reform", "2 Years Following Reform"),lty=c(2,1),lwd=c(1,1),cex=.6,bg="white")
dev.off()


################################## ################################## 
################# Figure B8 - Voter registration

pdf("Figures/FigureB8.pdf",width=8, height=4)
par(mfrow=c(1,2))
### Belgium
data <- read.dta("./Datasets/Other/Belgium_master.dta")
data<- subset(data,pctnoneunoncit >= quantile(data$pctnoneunoncit, c(.5),na.rm=T))  # Above median
data <- subset(data,data$year == 2005)
data <- subset(data,data$twoyear_growth > 0)
data <- subset(data,data$pop > 1000)  # registration rates fluctuate for small munis
data <- na.omit(data.frame(data$twoyear_growth,data$registered))
x <- (data$data.twoyear_growth)*100
y <- data$data.registered*100
general_loess(x,y,.99,"% Change in Per Head Spending, 2005 - 2003","Foreign Registration Rate (%)",c(10,30),"Belgium",.05,.95)

##### Switzerland
# No data available from Fribourg
data <- read.dta("./Datasets/Other/Swiss_master.dta")
data <- subset(data,(data$cantonid==22 & data$year == 2005) | (data$cantonid==25 & data$year==2007))
swissdata <- subset(swissdata,swissdata$pctforeign >= quantile(swissdata$pctforeign, c(.5),na.rm=T)) # Above median
data <- subset(data,data$twoyear_growth > 0)
data <- subset(data,data$pop > 1000)  # registration rates fluctuate for small munis
mdata <- na.omit(data.frame(data$twoyear_growth,data$foreign_participation))
x <- mdata$data.twoyear_growth*100
y <- mdata$data.foreign_participation
general_loess(x,y,.99,"% Change in Per Head Spending, 2 years after Reform","Foreign Registration Rate (%)",c(25,40),"Switzerland",.05,.95)
dev.off()


################################## ################################## 
################# Figure B9 - GAM 

pdf("Figures/FigureB9.pdf",width=8, height=8)

par(mfrow=c(3,3),mar=c(4,4,3,.5))

trimmed_gam <- function(data,title){
data <- subset(data,(S4_pctforeign > quantile(data$S4_pctforeign,.025,na.rm=T)) & (S4_pctforeign < quantile(data$S4_pctforeign,.975,na.rm=T)))
data <- subset(data,leftshare_begin < quantile(data$leftshare_begin,.975,na.rm=T)   & leftshare_begin > quantile(data$leftshare_begin,.025,na.rm=T) )
data <- subset(data,S4_log_social_head  < quantile(data$S4_log_social_head,.99,na.rm=T)   & S4_log_social_head > quantile(data$S4_log_social_head,.01,na.rm=T) )
data$S4_log_social_head  <- data$S4_log_social_head - mean(data$S4_log_social_head,na.rm=T)
mod.gam <- gam(S4_log_social_head ~ s(S4_pctforeign,leftshare_begin,sp=40),data=data)
vis.gam(mod.gam,plot.type="contour",color="topo",main=title,xlab="Change in % Foreign",ylab="Left Share",nCol=700,cex.lab=.8)
rug(data$S4_pctforeign)
rug(jitter(data$leftshare_begin,.01),side=2)
}

##### Belgium
datab <- read.dta("./Datasets/Crossnational/Belgium.dta")
data <- datab %>% group_by(id) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_cpas - lag(log_cpas,4))  %>% arrange(id,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Belgium")

##### Denmark
datab <- read.dta("./Datasets/Crossnational/Denmark.dta")
data <- datab %>% group_by(tempid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_net_social_head - lag(log_net_social_head,4)) %>% arrange(tempid,year)
data <- subset(data,year==2005 & !is.na(leftshare_begin) )
trimmed_gam(data,"Denmark")

##### Netherlands 
datab <- read.dta("./Datasets/Crossnational/Netherlands.dta")
data <- datab %>% group_by(tempid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_net_social_head - lag(log_net_social_head,4)) %>% arrange(tempid,year)
data <- subset(data,year==2009 & !is.na(leftshare_begin))
trimmed_gam(data,"Netherlands")

##### Norway 
datab <- read.dta("./Datasets/Crossnational/Norway.dta")
data <- datab %>% group_by(tempid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_net_social_head - lag(log_net_social_head,4)) %>% arrange(tempid,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Norway")

##### Sweden
datab <- read.dta("./Datasets/Crossnational/Sweden.dta")
data <- datab %>% group_by(tempid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_net_social_head - lag(log_net_social_head,4)) %>% arrange(tempid,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Sweden")

##### France
datab <- read.dta("./Datasets/Crossnational/France.dta")
data <- datab %>% group_by(tempid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_social_head - lag(log_social_head,4)) %>% rename(leftshare_begin=total_left) %>% arrange(tempid,year)
data <- subset(data,year==2009 & !is.na(leftshare_begin))
trimmed_gam(data,"France")

##### Italy
datab <- read.dta("./Datasets/Crossnational/Italy.dta")
data <- datab %>% group_by(id) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_social_head- lag(log_social_head,4)) %>% arrange(id,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Italy")

##### Spain
datab <- read.dta("./Datasets/Crossnational/Spain.dta")
data <- datab %>% group_by(mid) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_social_head - lag(log_social_head,4)) %>% arrange(mid,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Spain")

##### Switzerland
datab <- read.dta("./Datasets/Crossnational/Swiss_German.dta")
data <- datab %>% group_by(bfsnr) %>% mutate(S4_pctforeign = pctforeign - lag(pctforeign,4)) %>% mutate(S4_log_social_head = log_net_welfare_head - lag(log_net_welfare_head,4)) %>% arrange(bfsnr,year)
data <- subset(data,year==2008 & !is.na(leftshare_begin))
trimmed_gam(data,"Switzerland")
dev.off()